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Light harvesting hy photosystem II (PSII) in plants is highly efficient and accli¬ 
mates to rapid changes in the intensity of sunlight. However, the mechanisms 
of PSII light harvesting have remained experimentally inaccessible. Using a 
structure-hased model of excitation energy flow in 200 nanometer (nm) x 200 
nm patches of the grana membrane, where PSII is located, we accurately simu¬ 
lated chlorophyll fluorescence decay data with no free parameters. Excitation 
movement through the light harvesting antenna is diffusive, but becomes subd- 
iffusive in the presence of charge separation at reaction centers. The influence 
of membrane morphology on light harvesting efficiency is determined by the 
excitation diffusion length of 50 nm in the antenna. Our model provides the 
basis for understanding how nonphotochemical quenching mechanisms affect 
PSII light harvesting in grana membranes. 
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Main Text 


Photosynthetic light harvesting is a multi-scale process that spans from tens of femtoseconds 
to minutes, and from angstroms to hundreds of nanometers (1). In green plants, photosynthetic 
light harvesting is performed by photosystems I and II, which are located in separate regions 
of the thylakoid membrane (2). Photosystem II (PSII), which provides the electrons that drive 
photosynthesis, is located in the grana stacks of the thylakoid. PSII light harvesting starts when 
pigments bound to antenna proteins absorb sunlight. The nascent excitation energy is trans¬ 
ferred to reaction centers (RC) where it is converted to chemical energy via charge separation. 
PSII is both an efficient and adaptable light harvesting material. At maximum light harvesting 
capacity, the average time for excitation capture by PSII (~300 ps (3, 4)) is significantly less 
than the excited state lifetime of pigments in the grana (~2 ns (5, 6)), which results in a >80% 
efficiency for converting absorbed sunlight into chemical energy (7). In addition, PSII accli¬ 
mates to changes in sunlight intensity (5), wavelength (9), and downstream metabolism {10), 
which occur on the millisecond to minutes timescale. The final outcome of excitation in the 
grana is determined by the interplay between the rates of excitation energy transfer through 
the antenna and processes that irreversibly quench excitation (e.g. fluorescence, non-radiative 
processes, and productive photochemistry in the RC). Time-resolved chlorophyll fluorescence, 
which is currently the primary technique for measuring PSII light harvesting in a variety of en¬ 
vironmental conditions in vivo {11), does not have sufficient spatiotemporal resolution to fully 
characterize the light harvesting kinetics of PSII (72). Thus a model for PSII light harvesting 
that accurately captures its underlying dynamics is needed to understand the physical principles 
that give rise to its functionality. Such principles are essential for understanding the molecular 
mechanisms of photosynthesis and could be used to engineer artificial light harvesting systems 
with properties that mimic PSII {13-15). 
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Current models of light harvesting in PSII treat the rates of energy transfer and trapping as 
parameters whose values are determined by fitting to chlorophyll fluorescence data (4,16,17). 
Conceptually, these models mainly fit within the canonical “lake” and “puddle” models (18), 
which are currently used to describe the dynamics of energy transfer and trapping by PSII on the 
hundreds of picoseconds time scale in grana membranes. In the lake model, each excitation tra¬ 
verses the grana membrane sufficiently quickly to access all RCs equally. In the puddle model, 
excitation movement is spatially limited to the smallest photosynthetic unit, which is thought 
to be a photosystem II supercomplex (PSII-S) and a few surrounding major light harvesting 
complex II (LHCII) antenna (Fig. lA, inset) (19). PSII-S consists of two reaction centers and 
a small number of antenna proteins (20). Up to now, models with free parameters for energy 
transfer and trapping in PSII have been considered accurate if they can fit chlorophyll fluo¬ 
rescence decay data well. However, multiple models that implicitly assume either the lake or 
puddle model can fit chlorophyll fluorescence decay data equally well (4,16). As a result, it 
remains unclear whether the fitted rates found by these methods are physically meaningful and, 
more broadly, to what extent either the lake or puddle models are appropriate for describing 
PSII light harvesting. It is well-established that simulating the ultrafast excitation energy trans¬ 
fer dynamics in individual pigment-protein complexes requires a structure-based approach that 
correctly accounts for each chromophore (21,22). We previously showed that this approach is 
required even in a small group of pigment-protein complexes in our model of light harvesting 
within isolated PSII-S (23). Here, we demonstrate that an accurate model for PSII light har¬ 
vesting in grana membranes must be based on a correct description of the picosecond dynamics 
that occur at the nanometer scale of individual pigments. 

We constructed a parameter-less model for PSII light harvesting that is firmly based on an 
accurate physical description of the nanoscopic energy transfer dynamics in grana membranes. 
We performed Monte Carlo simulations on 200 nm x 200 nm patches of the grana membrane 
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(24) (see Methods, below) to generate examples of the mixed (Fig. lA) and segregated (Fig. IB) 
organizations previously observed (2, 26). In the segregated membrane, PSII-S and LHCII 
separate into so-called PSII-S crystalline arrays and LHCII aggregates. The Monte Carlo model 
contains a small number of energetic interactions based on in vivo phenomenology and provides 
plausible locations of LHCII in grana, which are difficult to visualize using existing imaging 
techniques (24). We superimposed the crystal structures of the chlorophyll pigments in PSII- 
S (20) and LHCII (27) on these simulations to establish the locations of all of the chlorophylls in 
a grana patch (Fig. lA, inset). In previous work on PSII-S, we grouped chlorophylls into highly 
energetically coupled clusters called domains and demonstrated that it is sufficient to describe 
the excitation dynamics at the domain level (23). We assume that the kinetics within LHCII 
and PSII-S are the same as calculated previously on isolated complexes (23). Inhomogeneously 
averaged rates of energy transfer between domains on different complexes were calculated using 
Generalized Forster theory (Methods). We modeled photochemistry in the reaction centers with 
two phenomenological “radical pair” (RP) states. The rate constants between the RP states are 
the same as those established in previous work on the PSII-S (23). Using our method we can 
routinely construct a rate matrix for grana membrane patches that contain more pigments (up to 
10^) than any natural light harvesting system previously modeled. 

Our model accurately simulates chlorophyll fluorescence data and demonstrates that ex¬ 
tracting the amplitude and lifetime components from a simple fit does not capture the complex 
kinetics of PSII light harvesting. The comparison of the simulation of the mixed membrane 
(black solid line) with experimental data from thylakoids in ref. (3) (red dotted line) is shown in 
Fig. 1C. We expected the mixed membrane to be the dominant morphology of the measured thy¬ 
lakoids based on the conditions in which the plants (Arabidopsis thaliana) were grown (3,28). 
The excellent agreement of the simulation with the experimental data, which was not guaranteed 
a priori, suggests that our parameter-free model correctly describes excitation energy transport 
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during PSII light harvesting. The simulated deeay ean be ht well to a sum of a few exponentials 
(Fig. 1C, inset, green bars), as is frequently done to extraet the amplitude and lifetime eompo- 
nents {11). However, the ht does not eapture the eomplex distribution ealeulated from the rate 
matrix (Fig. 1C, inset, blaek bars). Fit-based models that are only sufheiently eomplex to ht 
ehlorophyll fluoreseenee deeays well will not aeeurately eharaeterize the underlying PSII light 
harvesting dynamies. 

Our simulation of the spatiotemporal dynamies of ehlorophyll exeitation in the grana al¬ 
lowed us to examine the assumptions underpinning the lake and puddle models. Upon uniform 
initial exeitation of either the mixed or segregated membranes, the lake model prediets that the 
exeitation distribution will quiekly reaeh a steady-state spatial distribution. However, we did 
not observe a steady-state distribution for either the mixed (Movie SI) or segregated (Movie 
S2) membranes. The puddle model prediets that there are a few (2-3) different types of trapping 
dynamies at reaetion eenters throughout the membrane. On the eontrary, we observed trapping 
dynamies (Movies S3 and S4) that vary eonsiderably between different PSITS. Neither the lake 
nor the puddle models aeeurately deseribe exeitation dynamies on the hundreds of pieoseeonds 
timeseale and hundreds of nanometers length seale of the grana membrane. 

To understand how exeitation moves through the grana, we simulated exeitation energy flow 
from single pigment-protein eomplexes in LHCII aggregates, PSII erystalline arrays, and mixed 
membranes (Movies S5-S9). To determine the effeet of eharge separation on exeitation move¬ 
ment, we simulated the PSII erystalline arrays and mixed membranes both with and without 
the radieal pair (RP) states in the reaetion eenters. We quantified exeitation transport in these 
simulations by ealeulating the time-dependenee of the varianee of the exeitation probability 
distribution (Fig. 2A). Transport aeross the grana was well deseribed by fitting the equation 

(J^{t) — cr^(O) = (1) 
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A Mixed membrane B Segregated membrane 



Figure 1: Structure-based modeling of energy transfer is required to accurately simu¬ 
late PSII light harvesting. (A) and (B) The representative mixed and segregated membrane 
morphologies generated using Monte Carlo simulations and used throughout this work. PSII 
supereomplexes (PSII-S) are indieated by the teal diseoreetangles, while major light harvesting 
antenna eomplexes (LHCII) are indieated by the green eireles. The segregated membrane forms 
PSITS erystalline arrays and LHCII aggregates. As shown sohematieally in the inset in (A), ex¬ 
isting crystal structures of PSII-S (20) and LHCII (27) were overlaid on these membrane patches 
to establish the locations of all chlorophyll pigments. The teal and green dashed lines outline the 
excluded area associated with PSII-S and LHCII in the Monte Carlo simulations, respectively. 
The chlorophyll pigments are indicated in blue, while the protein is depicted by the cartoon 
ribbon. PSII-S is a 2-fold symmetric dimer of pigment-protein complexes that are outlined by 
black lines. LHCII-s, CP26, CP29, CP43, and CP47 are antenna proteins, while RC indicates 
the reaction center. The inhomogeneously-averaged rates of energy transfer between strongly- 
coupled clusters of pigments were calculated using Generalized Forster theory. (C) Simulated 
fluorescence decay of the mixed membrane (solid black line) and the PSII-component of exper¬ 
imental fluorescence decay data from thylakoid membranes from ref. (3) (red, dotted line). The 
inset shows the lifetime components and amplitudes of the simulated decay as calculated using 
our model (black bars) or by fitting to three exponential decays (green bars). 
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where cT^(t) is the variance at time t, cr^(O) is the variance of the initial distribution of excitation, 
and A and a are fit parameters (Fig. S4B). For a = 1, transport is diffusive, while subdiffusive 
transport results if a is significantly less than 1. For the LHCII aggregate (Fig. 2A, green tri¬ 
angles) and for the mixed membranes and PSII crystalline arrays without RP states (Fig. 2A, 
red and blue dashed lines, respectively), a ~ 1 and thus transport within the antenna can be 
considered diffusive. However, transport for both the mixed and PSII crystal cases with RP 
states (Fig. 2A, red and blue solid lines, respectively) is sub-diffusive {a < 0.75). Subdiffu- 
sivity can occur when the energetic differences between sites is on the order of or greater than 
ksT. We calculated the AG for the RC—>RP1 (radical pair 1 ) step to be -SAksT on the basis 
of our previously published rates (23). Thus, RPl serves as an energetic trap that causes sub- 
diffusive transport. The energy transfer rates in our model are averaged over inhomogeneous 
realizations, which could mitigate the slow-down of diffusion that occurs when the width of the 
inhomogeneous distribution is greater than ksT (29). The largest standard deviation of exciton 
energies across an inhomogeneous distribution in our model is 107 cm“\ which is significantly 
less than the value of ksT at room temperature (210 cm“^). We note that the diffusive picture of 
excitation energy transport elaborated here is consistent with predictions made from the PSITS 
model (23) and suggests that, in the grana, excitation energy flows neither “directionally” nor 
energetically downhill on “preferred pathways” (30). 

We calculated the excitation diffusion length, L^, and the diffusivity, D, to reduce the 
dynamics observed above to single parameters that can be used to qualitatively understand light 
harvesting behavior in grana. 

a\t) - a\0) = 4Dt = L'^ (2) 

A plot of D as a function of time is shown in Fig. S4A. In the case of subdiffusive transport, 
the diffusivities decrease in time as excitations are held by the RPl trap state. The diffusion 
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Legend 

— Mixed with RP 

— PSii-S c.a. with RP 

— Mixed no RP 

-♦-PSii-S c.a. no RP 



Figure 2: Excitation transport in grana membranes. Simulation of excitation movement 
in the five grana membrane configurations shown in the legend: mixed membrane with and 
without the radical pair (RP) states in the reaction center, PSII-S crystalline array (PSII-S c.a.) 
with and without the RP states, and LHCII aggregate (LHCII agg). In each case, excitation 
was initiated on a single pigment-protein complex. (A) The change in the spread of excitation 
over time. The diffusion exponent a (see text for details) is shown on the right of the plot. 
(B) Fraction of surviving excitation as a function of net displacement L (eq. 2) from the initial 
starting point. The dashed line, where the fraction of surviving excitation is 1/e, demarcates the 
excitation diffusion length (Ljj). The dimensions of some of the configurations were too small 
to calculate an L^), so linear extrapolation was used to approximate it (line segments that do not 
include markers). 
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constants for the three cases ranged from 1—5 x 10“^ cm^/sec. These values agree with previous 
experimental work, which used singlet-singlet annihilation measurements of chloroplasts to 
suggest a lower bound of D ^ 10“^ cm“^/sec (31). The fraction of excitation remaining as a 
function of the net spatial displacement, L, is shown in Fig. 2B. is defined, by convention, as 
the minimum net displacement in one dimension achieved by 37% of the excitation population. 
The Lf) for the mixed membrane and the PSII crystal with radical pair states were ~25 nm 
and ~15 nm, respectively, though there was some variation depending on the starting point of 
excitation (Fig. S5). The for the LHCII aggregate, mixed membrane, and PSII crystalline 
array without radical pair states was ~50 nm. The values of and D in the antenna compare 
favorably with the values observed thus far in organic thin films (32) and quantum dot arrays 
(29). 

Using the excitation diffusion length (Lc), we explored the longstanding question of how 
membrane morphology influences the efficiency with which PSII converts absorbed light into 
chemical energy, the photochemical yield (2, 30). We calculated the photochemical yield in 
both the mixed and segregated cases to be 0.82 and 0.70, respectively, upon spatially uniform 
excitation across the membranes. The calculated value of 0.82 for the mixed membrane is 
in excellent agreement with the estimate of 0.83 derived from chlorophyll fluorescence yield 
measurements (7). Previous work suggested that the reduced light harvesting efficiency of seg¬ 
regated membranes was due to an inability of excitation initiated in LHCII aggregates to drive 
photochemistry (3), resulting in the “disconnection” of LHCIIs from reaction centers. To ex¬ 
plore this hypothesis we initiated excitation on each LHCII in both the mixed and segregated 
membranes and calculated the resulting photochemical yield, $lhcii (Fig- 3A). Both mem¬ 
branes have a wide distribution of $lhcii (Fig. 3B). (<Flhcii) for the segregated membrane, 
0.49, is significantly lower than that for the mixed membrane, 0.75. This difference can be 
explained most simply by the relative diffusion lengths in each membrane environment. In the 
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Figure 3: Influence of grana membrane morphology on photochemical yield. (A) Excitation 
was initiated at each LHCII in both the mixed (left) and segregated (right) morphologies. The 
color of the LHCII indicates the fraction of excitation that results in productive photochemistry 
(‘^’lhciu see colorbar on far right) as simulated with our model. (B) Histograms representing 
the distribution of <hLHCii for the mixed and segregated membranes using the coloration from 
(A). 


mixed membrane case, nearly all LHCIIs exist within an (25 nm) of a PSII-S. In the seg¬ 
regated membrane, however, the LHCII aggregate(s) have a diameter comparable to (50 
nm), resulting in substantial loss of excitation prior to reaching a PSII-S. Clearly, LHCIIs in 
the segregated membrane show reduced light harvesting function; however, the LHCIIs do not 
completely disconnect (<I>lhcii ~ 0). To disconnect an LHCII aggregate from the neighboring 
PSII-S requires the aggregate diameter to exceed 100 nm, which is unlikely to occur given that 
a grana membrane disc has a diameter of ~400 nm (2). The physiological benefit of segre¬ 
gation remains unclear and may be explained by other aspects of photosynthesis, as has been 
proposed (33,34). 

Our model accurately describes the energy transfer network of the grana membrane and 
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thus paves the way for understanding PSII light harvesting in all environmental conditions. In 
intense sunlight, the rate of light absorption by the antenna exceeds the rate of trapping by 
reaction centers, and the need for photoprotection arises. The mechanisms underlying the non¬ 
photochemical quenching of excess excitation in the antenna remain controversial (5). Most 
of the hypothesized mechanisms of quenching are based, by necessity, on measurements per¬ 
formed on isolated pigment-protein complexes {11). However, it has been unclear whether the 
mechanisms observed in vitro occur in vivo. Our model provides a unified framework for un¬ 
derstanding to what extent picosecond dynamics observed in vitro explain in vivo chlorophyll 
fluorescence data. More broadly, PSII is connected to the rest of photosynthesis through the 
electron transport chain and the pH gradient across the thylakoid membrane (35,36). Our model 
could be extended to include the unappressed portion of thylakoid membranes, which contains 
photosystem I, and integrated into models of the entire thylakoid to fully model photosynthesis 
in plants. 

Methods 

Overview 

There has been an extensive theoretical development of energy transfer models in photosyn¬ 
thetic complexes composed of a few proteins. Excitation dynamics can be calculated directly 
from the Hamiltonian of a pigment-protein complex (37). However, exact calculations of the dy¬ 
namics using the Hamiltonian are computationally costly, and, while feasible on complexes with 
~100 pigments, are impractical for simulating dynamics in the photosystem Il-containing por¬ 
tions of the thylakoid membrane, which contain hundreds of complexes and >10,000 pigments. 
Another common technique for calculating excitation dynamics is by perturbative treatments of 
the Hamiltonian, such as the Generalized Forster and Redfield methods (21). In photosystem II 
(PSII), the energetic coupling between pigments spans a wide range, such that neither method is 
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appropriate on its own. Modified Redfield theory can interpolate between the strong and weak 
coupling regimes, but it does not account for dynamic localization in which coupling between 
the phonon modes of the protein and the pigments’ excited states restrict the delocalization of 
excitonic states. Therefore, energy transfer in PSII has been treated using a combination of the 
Modified Redfield and Generalized Forster theories (23,38,39). Transfer within tightly coupled 
clusters of pigments is treated by Modified Redfield theory and transfer between these clusters 
is treated using Generalized Forster theory. 

In previous work on PSII supercomplexes (PSITS), we used the Modified Redfield/Generalized 
Forster (MR/GF) method in combination with a novel approach for defining highly coupled 
clusters of chlorophylls, or domains, that optimizes the separation of timescales between intra- 
and inter-domain transport (23). Using these domain definitions we were able to coarse-grain 
the excitation dynamics at the domain level. Recent work using a more exact method (ZOFE) 
for calculating energy transfer dynamics has further validated our MR/GF simulations (40). 
Another group has performed HEOM calculations on a quadrant of the PSITS (41). A compar¬ 
ison of our MR/GE simulations with these calculations gives good agreement when appropriate 
domain definitions are used. 

We have used the same domain definitions in our simulations of energy transfer in thy- 
lakoid membranes as were used for modeling PSITS. We assumed that the timescales of en¬ 
ergy transfer between domains on different complexes in the membrane are slow relative to 
the intra-domain timescales (<1 ps“^) calculated previously (23) and used Generali z ed Edrster 
theory to calculate the rates between domains on different complexes (see Excitation energy 
transfer section below). The locations of all of the PSITS and major light harvesting complex II 
(EHCII) antenna in a grana membrane were calculated on the basis of Monte Carlo simulations 
of thylakoid membranes (24) (see Monte Carlo simulations of paired grana membranes section 
below). The crystal structures of PSITS (20) and EHCII (27) were overlaid on top of these 
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simulations (see Chlorophyll configurations section below). 

The dynamics of PSII light harvesting can be represented by a kinetic network composed of 
Ist-order rate constants. The master equation formalism can be used to calculate the dynamics 
of excitation: 


P{t)=KP{t), (3) 

where K is a rate matrix containing the first order rate constants of excitation transfer between 
all compartments (which include all domains) in the network, and P{t) is the vector of compart¬ 
ment populations (11,42). Within this framework, we use simple models for charge separation 
and trapping when excitation reaches the reaction center (see Electron transfer model section 
below). While there are more complex models for electron transfer within the reaction cen¬ 
ter (43), the two-compartment model used here was parameterized on PSII-S with different 
antenna sizes (23) and is sufficient to determine the overall timescales of this process. The 
non-radiative rate constant of decay from each domain was (2 ns)“^ (23). The fluorescence rate 
constant for each exciton was scaled by its transition dipole moment squared with the average 
fluorescence rate constant across all excitons set to (16 ns)“^ (23). 

Solving eq. 1 for P(t), which contains the dynamics of excitation population on all com¬ 
partments in the network, gives 


P(t) = CP^C-^P(O), (4) 

where C is a matrix which contains the eigenvectors of iT, A is a diagonal matrix containing 
the eigenvalues of K, and P(0) is the initial vector of populations. Using eq. 2, we calculate 
chlorophyll fluorescence decays, the dynamics of excitation over time, and the yields of the 
different dissipation processes available to excitation - non-radiative decay, fluorescence, and 
productive photochemistry in the reaction centers (see Calculations of P(t) section below). 
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Monte Carlo simulations of paired grana membranes 

Grana-scale pigment-protein eomplex eonfigurations were generated via eomputer simulations 
of an extension of the model presented in ref. (24). Briefly, dise-shaped partieles L representing 
LHCII complexes and rod-shaped particles P representing so-called C 2 S 2 PSII-S in 2d layers a 
and (3 interacted via hard-core repulsive interactions, plus the attractive energetic potentials 



(5) 



( 6 ) 


with Cstack = 4 k-^T, Cm = 2 k-^T, Aagg = 1.15, and all other parameters as in ref. (24). 

The potentials in Eqs. and are motivated in ref. (24). The square-well attraction in 
Equationj^acts a phenomenological, non-specific energetic driving force for EHCII aggregation 
(44). In the “mixed” condition, Cagg = 0 k-^T, while in the “segregated” condition, Cagg = 1 
ksT. 

Canonical ensemble Metropolis Monte Carlo simulations of this model were performed as 
in ref. (24). A ratio of free EHCII to PSITS particles of 6 was chosen to match the conditions 
of ref. (3), and a particle packing fraction of 0.75 was chosen to match typical grana conditions. 
Thus, 64 PSITS particles and 384 free EHCII particles were initialized in each of two square 
boxes of side length 200 nm. Simulations were equilibrated with periodic boundary conditions 
for at least 15M Monte Carlo sweep steps. 

Chlorophyll configurations 

Representative configurations from the stroma-side-up layers of the Monte Carlo simulations 
were selected for analysis. Eor each configuration, chlorophyll coordinates were assigned for 


14 



each pigment-protein particle by aligning the center and axis of rotation of the chlorophyll coor¬ 
dinates from refs. (20,27) to the center and axis of rotation of the simulated particle (Fig. S3). As 
discussed and motivated in ref. (23), we have substituted the structure of an LHCII monomer in 
the place of the minor light harvesting complexes in PSII-S. Because simulated LHCII particles 
were radially symmetric, an axis of rotation was randomly selected for each LHCII particle. 

In a small number of cases across the membrane, pigments belonging to one complex enter 
into the excluded area of a different complex. This overlap originates from a mismatch between 
the idealized geometries used for the Monte Carlo simulations and the real pigment structure 
as shown in Fig. S2A. As a result of protein overlap, there exist a small number of analogously 
high transferrates within the membrane rate matrix. In order to understand the influence of such 
rates on the overall description of transport in the membrane, we have removed transfer rates 
exceeding certain thresholds from the rate matrix and re-calculated the resulting fluorescence 
decay. As can be seen in Fig. S2B, this does not alter the fluorescence decay dynamics. 

Excitation Energy Transfer Model 

The rate of energy transfer from a donor domain d to an acceptor domain a, (eq. 3), is 
typically calculated for each inhomogeneous realization of a pigment-protein complex using 
generalized Forster theory (23,38,39). kf^^ is the Boltzman-weighted (eq. 4) sum of the rates 
from the excitons in d (\M) G d) to the excitons in a (|A^) G a). The rates between excitons 
are calculated using the generalized Forster equation (eq. 5-6), where dtAM(t)F^(t) is 
the overlap integral, \Vm,n\'^ is the electronic coupling between the two excitons, and 
is the coefficient of the Mth exciton on the site jj,. Calculation of some observable of the 
system, such as the fluorescence lifetime, is usually done for each realization and then averaged. 
However, the thylakoid membrane contains >10,000 chlorophylls. Generating hundreds of 
realizations of the population coefficients and the overlap integrals for the membrane is not 
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only computationally intensive, but may not necessary to accurately describe the dynamics of 
the system at the protein length scale. 
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dtAM{t)F^{t) 


( 10 ) 


\Vm,n\^ = (11) 

Using the inhomogeneously averaged rate matrix gives the correct overall dynamics and is 
representative of the energetics of PSII, while greatly increasing the computational efficiency. 
It is computationally feasible to calculate a single rate matrix for the intact thylakoid mem¬ 
brane using MR/GF theory through the use of supercomputers and patience. This calculation, 
however, is likely to be unnecessary, because a much more computationally efficient method 
is available for calculating the inhomogeneously averaged rate matrix, which gives the correct 
overall dynamics, as shown by the fluorescence lifetime comparison in Fig. SI for a PSITS. 
While the dynamics appear to be accurate, this averaging flattens out the energy landscape and 
will speed up diffusion, as has been shown in quantum dot arrays (29). The standard devi¬ 
ation of the exciton energies across the ensemble of inhomogeneous realizations in PSII are 
significantly less than ksT, suggesting that excitation movement through the grana membrane 
should not be greatly affected by averaging. Nonetheless, ongoing developments in making 
exact calculations of energy transfer more scalable (45) will hopefully enable such calculations 
on membranes in the future. 
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Directly calculating the inhomogeneous average transfer rates between domains substan¬ 
tially reduces the computational burden because a single PSII-S contains all of the different 
domains that occur throughout the grana membrane. The most computationally demanding 
components for determining the rate matrix for a given inhomogeneous realization are the over¬ 
lap integral in eq. 5 and the matrix to transform from the site to the exciton basis in eq. 6. We 
have assumed that domain definitions, overlap integrals, and the transformation matrices for 
PSITS and LHCII are the same at the membrane level as they are in isolated complexes. As 
a result these terms have been computed previously for several hundred realizations of inho¬ 
mogeneous broadening for the four different types of transfers (from LHCII to LHCII, PSII to 
PSII, LHCII to PSII and PSII to LHCII) in the membrane (23). We tabulated these values and 
combined them with a calculation of the electrostatic coupling between sites (the only term that 
must be calculated for each domain-to-domain rate in the membrane). By using our tabulated 
values for the population matrices and overlap integrals, we have reduced the computational 
time by 3-4 orders of magnitude, turning an otherwise burdensome calculation into one that can 
be performed in less than 1 day on a single CPU. We note that this simplification works only 
because we calculate the inhomogeneously averaged rate between each domain. This allows for 
a much smaller sampling of the possible combinations of site energies then would be required 
for even a single inhomogeneous realization of the entire thylakoid membrane. 

Our algorithm for calculating {kf^^)iahom. real, for each domain pair in the membrane was 
as follows. We first checked if the two domains were within 60 A of each other. If yes, we 
proceeded to calculate the average rate; if no, the rate was set equal to 0. We used 60 A as a 
cutoff based on the distance dependence of rates in the PSII supercomplex. If the two domains 
were from the same complex, we used the inhomogeneously averaged rate of transfer from 
already performed generalized Forster/modified Redfield calculations (23). If not, we noted 
in which type of complex the donor domain and the acceptor domain were in, either LHCII 
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or PSII-S. We calculated using eq. 3-6 with pre-tabulated values of the overlap integral, 
the population matrices, and the energies of the excitons for a few hundred inhomogeneous 
realizations. The electronic coupling was uniquely calculated for that pair of domains 
using the ideal dipole approximation. We then averaged over these tabulated inhomogeneous 
realizations to get (fc^^)inhom. real.- 

Electron transfer model 

Both the identity of the primary donor and the kinetics and mechanism of charge separation in 
the PS II reaction center remain controversial {46). The lack of agreement between the various 
experimental results and theory has resulted in a variety of phenomenological and conceptual 
models being used to interpret experimental data {17,47). In previous work on PSII supercom¬ 
plexes, we used the simplest kinetic model that describes the processes known to occur in the 
PSII reaction center (23), 


( 12 ) 

Here, RC is the reaction center domain composed of the 6 pigments of the reaction center. The 
“radical pair” states RPl and RP2 and the rate constants knc, kcs, and kirr are used to model 
the electron transfer steps in the reaction center. RPl and RP2 are non-emissive states that do 
not have a direct physical analog with charge-separated states in the reaction center. Rather, 
this approach allows us to describe the overall process of a reversible charge separation step 
followed by an irreversible step and establish the approximate timescales of these events relative 
to energy transfer in the light harvesting antenna. The rates were previously parameterized 
using fluorescence decays of PSII supercomplexes of different sizes {19), with tcs = 0.64 ps. 
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'Trc = 160 ps, and Tj^r = 520 ps (/c = 1/r) (23). 

Calculations of P{t) 

Using the definitions established in the previous seetions, we ean write K as follows: 


udom 

i,j < 7^ j 

(13) 

k 

i e RC,j e RPl,j > N^om 

(14) 

k 

i e RPl,j e RC,i > A^dom 

(15) 

k- 

i e RPl,j G RP2, i,j > Ndom 

(16) 

^dump? 

'1 — P^dora ) j -^comp 1 

(17) 

UI 

i ^ dom) j ~ A^comp 

(18) 

\ ^ udom 

i = j. 

(19) 




where Kij is a matrix element of row i and eolumn j, iVdom is the total number of domains, and 
IV^comp is the total number of eompartments (iVdom plus all RPl, RP2, FI, and dump eompart- 
ments). The sizes of the K for the mixed and segregated membranes modeled (Fig. 1A-B, main 
text) were both 10882 x 10882, which meant that calculating the eigenvalues and eigenvectors 
required for simulating the dynamics of excitation (eq. 2) required the use of supercomputers 
with ~30 GB of memory. 

Fluorescence decays were calculated using the equations described in ref. (23). P(0) in all 
cases was that for ChlA excitation. The photochemical yield was calculated by summing over 
the populations in all RP2 states at t = 1 sec. 
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